

```r
#INTRO: LOADING PACKAGES ETC----
library(foreign)
library(reshape2)
library(data.table)
```

```
## 
## Attaching package: 'data.table'
```

```
## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt
```

```r
library(stargazer)
```

```
## 
## Please cite as:
```

```
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
```

```
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
```

```r
library(MatchIt)
#library(nonrandom)
library(ggplot2)
#library(gdata1)
library(gmodels)
library(gridExtra)
library(haven)
library(plyr)
#library(car)
#library(xlsx)
library(tidyr)
```

```
## 
## Attaching package: 'tidyr'
```

```
## The following object is masked from 'package:reshape2':
## 
##     smiths
```

```r
library(reshape)
```

```
## 
## Attaching package: 'reshape'
```

```
## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths
```

```
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
```

```
## The following object is masked from 'package:data.table':
## 
##     melt
```

```
## The following objects are masked from 'package:reshape2':
## 
##     colsplit, melt, recast
```

```r
#library(countrycode)
#library(dGlyr)
library(readr)
library(gnm)
```

```
## Error in library(gnm): there is no package called 'gnm'
```

```r
library(MNP)
```

```
## Error in library(MNP): there is no package called 'MNP'
```

```r
library(nls2)
```

```
## Error in library(nls2): there is no package called 'nls2'
```

```r
library(nlstools)
```

```
## Error in library(nlstools): there is no package called 'nlstools'
```

```r
library(stats)
library(dplyr)
```

```
## 
## Attaching package: 'dplyr'
```

```
## The following object is masked from 'package:reshape':
## 
##     rename
```

```
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
```

```
## The following object is masked from 'package:gridExtra':
## 
##     combine
```

```
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
```

```
## The following objects are masked from 'package:stats':
## 
##     filter, lag
```

```
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
```

```r
library(matrixStats)
```

```
## 
## Attaching package: 'matrixStats'
```

```
## The following object is masked from 'package:dplyr':
## 
##     count
```

```
## The following object is masked from 'package:plyr':
## 
##     count
```

```r
library(miceadds)
```

```
## Loading required package: mice
```

```
## 
## Attaching package: 'mice'
```

```
## The following object is masked from 'package:stats':
## 
##     filter
```

```
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
```

```
## * miceadds 3.13-12 (2022-05-30 15:14:07)
```

```r
library(minpack.lm)
library(lfe)
```

```
## Error in library(lfe): there is no package called 'lfe'
```

```r
library(lubridate)
```

```
## 
## Attaching package: 'lubridate'
```

```
## The following object is masked from 'package:reshape':
## 
##     stamp
```

```
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
```

```
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
```

```r
library(AER)
```

```
## Error in library(AER): there is no package called 'AER'
```

```r
library(reshape)
library(reshape2)
library(tictoc)
library(xtable)
library(Formula)
library(mlogit)
```

```
## Loading required package: dfidx
```

```
## 
## Attaching package: 'dfidx'
```

```
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate
```

```
## The following object is masked from 'package:stats':
## 
##     filter
```

```r
library(plm)
```

```
## Error in library(plm): there is no package called 'plm'
```

```r
library(survival)
library(plotly)
```

```
## 
## Attaching package: 'plotly'
```

```
## The following object is masked from 'package:reshape':
## 
##     rename
```

```
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise
```

```
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
```

```
## The following object is masked from 'package:stats':
## 
##     filter
```

```
## The following object is masked from 'package:graphics':
## 
##     layout
```

```r
library(alpaca)
```

```
## Error in library(alpaca): there is no package called 'alpaca'
```

```r
library(lmtest)
```

```
## Loading required package: zoo
```

```
## 
## Attaching package: 'zoo'
```

```
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
```

```r
library(sandwich)
# library(mnlogit)

library(foreach)#parallel foreach
library(iterators)
library(parallel)
library(doParallel)

library(msm)
library(tikzDevice)
```

```
## Error in library(tikzDevice): there is no package called 'tikzDevice'
```

```r
library(ggpubr)#combine ggplot graphs together
```

```
## 
## Attaching package: 'ggpubr'
```

```
## The following object is masked from 'package:plyr':
## 
##     mutate
```

```r
library(estimatr)#lm_robust
```

```
## Error in library(estimatr): there is no package called 'estimatr'
```

```r
library(clubSandwich)#coef_test
```

```
## Error in library(clubSandwich): there is no package called 'clubSandwich'
```

```r
library(sjlabelled)#labelled variables (set_labels)
```

```
## Error in library(sjlabelled): there is no package called 'sjlabelled'
```

```r
library(openxlsx)#read.xlsx

library(labelled)#change labels to columns

library(texreg)#extract for lm.cluster
```

```
## Error in library(texreg): there is no package called 'texreg'
```

```r
options(width=2000)#set max characters shown
options(max.print=200000)#set max rows shown console

# !diagnostics suppress=Dsvr,Dgdm

rm(list=ls())
cat("\014")
```



```r
#READ DATA:----
setwd("~/Dropbox/Princeton/P29 reassessing leaders effect/_data")
```

```
## Error in setwd("~/Dropbox/Princeton/P29 reassessing leaders effect/_data"): cannot change working directory
```

```r
load("d29.RData")


#ASSIGN WEIGHTS:----

#assign Weight Country:
#calculate country (weighted by country population) weights:
# de = de %>% group_by(Es) %>% mutate(Resp = n()/6 )#number of respondent by Es
de$Pop=NA#country population in millions
de$Pop[de$Ec=="AUS"]=25
de$Pop[de$Ec=="AUT"]=9
de$Pop[de$Ec=="CAN"]=37
de$Pop[de$Ec=="DEU"]=68#West Germany (84 All Germany)
de$Pop[de$Ec=="DNK"]=6
de$Pop[de$Ec=="ESP"]=47
de$Pop[de$Ec=="FIN"]=6
de$Pop[de$Ec=="GBR"]=61
de$Pop[de$Ec=="GRC"]=10
de$Pop[de$Ec=="IRL"]=5
de$Pop[de$Ec=="ISL"]=1
de$Pop[de$Ec=="ISR"]=9
de$Pop[de$Ec=="ITA"]=61
de$Pop[de$Ec=="NLD"]=17
de$Pop[de$Ec=="NOR"]=5
de$Pop[de$Ec=="NZL"]=5
de$Pop[de$Ec=="PRT"]=10
de$Pop[de$Ec=="SWE"]=10
de$Rwp=de$Pop#Respondent Weight Country
de$Rwp=de$Rwp/mean(de$Rwp)#set mean weight == 1

#calculate country weights (weights s.t. all elections same number of respondents):
de = de %>% group_by(Es) %>% mutate(Resp = n()/6 )#number of respondent by Es
de$Rwc=1000/de$Resp#Respondent Weight Country
de$Rwc=de$Rwc/mean(de$Rwc)#set mean weight == 1
summary(de$Rwc)
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3826  0.7022  0.8822  1.0000  1.2445  3.4331
```

```r
#set demographic weights to average 1 by election:
# (since I deleted respondents w/o Va,Pl,Ll, they unbalance the weights for that election
# in a way that is election specific, in this way I reset sum weghts to be 1 by election)
de = de %>% group_by(Es) %>% mutate(Rwd=Rwd/mean(Rwd))
summary(de$Rwd)
```

```
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.02013  0.98309  1.00000  1.00000  1.00000 11.38005
```

```r
#calculate country demographic weights:
# NB: this is equivalent to CSES Dataset Weight:
# | The derivative "Dataset Weight" (D1014) has been created so
# | that each election study in the dataset will contribute
# | equally to analyses of respondents, regardless of the number
# | of interviews in each election study.
de$Rwcd=de$Rwc*de$Rwd#Respondent Weight Country and Demographic
de$Rwcd=de$Rwcd/mean(de$Rwcd,na.rm=T)#set mean weight == 1
summary(de$Rwcd)
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0159  0.6837  0.8777  1.0000  1.1472 18.3894
```

```r
max(de$Rwcd)/min(de$Rwcd)
```

```
## [1] 1156.515
```

```r
de$Rwcp=de$Rwc*de$Rwp
de$Rwcp=de$Rwcp/mean(de$Rwcp,na.rm=T)#set mean weight == 1

de$Rwcdp=de$Rwc*de$Rwd*de$Rwp
de$Rwcdp=de$Rwcdp/mean(de$Rwcdp,na.rm=T)#set mean weight == 1

de=de[,-which(names(de) %in% c("Resp","Pop"))]#remove working columns


#CALCULATE Pl,Ll BY ELECTION:----
reg=de %>% group_by(Es) %>% summarise(mod = list(summary(clogit(Va ~ Pl+Ll +strata(Esalt), robust=T, weights=Rwcd, method="efron"))$coefficients[,1]))#create Pl and Ll by Es
length(unlist(reg))#366
```

```
## [1] 366
```

```r
Es=unlist(reg)[1:122]
Pl=unlist(reg)[123:366]
Pl=Pl[seq(1,244,2)]
Ll=unlist(reg)[123:366]
Ll=Ll[seq(2,244,2)]
ds=data.frame(Es,Pl,Ll)
ds$Pl=as.numeric(as.character(ds$Pl))
ds$Ll=as.numeric(as.character(ds$Ll))
ds$PlLl=ds$Pl-ds$Ll
ds$Ey=substr(ds$Es,5,8)
ds$Ey=as.numeric(ds$Ey)
ds$Ec=substr(ds$Es,1,3)#
ds$time=ds$Ey-1960
rm(reg,Es,Pl,Ll)

#merge Ed:
dss=de[,c(3,4)]
dss=dss %>% group_by(Es) %>% filter(row_number (Es) == 1)
ds=merge(ds,dss,by="Es")
ds$Ed=as.factor(ds$Ed)#(for graphs)

#merge NP:
dss=de[,c(3,7)]
dss=dss %>% group_by(Es) %>% filter(row_number (Es) == 1)
ds=merge(ds,dss,by="Es")
ds$NP=as.factor(ds$NP)#(for graphs)


#BASE FOR COUNTRY SLOPE FE:----
#Ec:Pl FE:
for (i1 in 1:length(sort(unique(de$Ec)))) {#creating dummy variables Ec:Pl (for FE)
  de[[paste0("Pl_",sort(unique(de$Ec))[i1])]]=ifelse(de$Ec==sort(unique(de$Ec))[i1],de$Pl,0)#de$Pl
}
#Ec:Ll FE:
for (i1 in 1:length(sort(unique(de$Ec)))) {#creating dummy variables Ec:Ll (for FE)
  de[[paste0("Ll_",sort(unique(de$Ec))[i1])]]=ifelse(de$Ec==sort(unique(de$Ec))[i1],de$Ll,0)#de$Ll
}
#create text for regression function:
txt=""
for (i1 in 1:length(sort(unique(de$Ec)))) {#creating function text for dummy variables Es:alt except for alt=1
  txt=paste(txt,paste0("Pl_",sort(unique(de$Ec))[i1]),sep="+")
  txt=paste(txt,paste0("Ll_",sort(unique(de$Ec))[i1]),sep="+")
}


#ALTERNATIVE SPECIFICATIONS SPLINE TWO PERIODS:----

#optimization:
cores=detectCores()
cl=makeCluster(cores[1])
registerDoParallel(cl)#set parallel computing clusters
TRses=foreach (t = 1:59, .packages='survival', .combine='c') %dopar% {#getting threshold (via max LogLik):
  de$time1=pmin(de$time,t)
  de$time2=pmax(t,de$time)
  CLC_ses_=clogit(Va ~ Pl:time1+Ll:time1+Pl:time2+Ll:time2
                     +Pl_AUS+Ll_AUS+Pl_AUT+Ll_AUT+Pl_CAN+Ll_CAN+Pl_DEU+Ll_DEU+Pl_DNK+Ll_DNK+Pl_ESP+Ll_ESP+Pl_FIN+Ll_FIN+Pl_GBR+Ll_GBR+Pl_GRC+Ll_GRC+Pl_IRL+Ll_IRL+Pl_ISL+Ll_ISL+Pl_ISR+Ll_ISR+Pl_ITA+Ll_ITA+Pl_NLD+Ll_NLD+Pl_NOR+Ll_NOR+Pl_NZL+Ll_NZL+Pl_PRT+Ll_PRT+Pl_SWE+Ll_SWE
                     +strata(Esalt), data=de, method="efron")#regression
  CLC_ses_$loglik[2]
}
stopCluster(cl)#stop parallel computing clusters
tr=which.max(TRses)#38
trCLC_ses=tr#threshold
trCLC_ses
```

```
## [1] 38
```

```r
#ALTERNATIVE SPECIFICATIONS SPLINE THREE PERIODS:----

#optimization:
tic()
cores=detectCores()
cl=makeCluster(cores[1])
registerDoParallel(cl)#set parallel computing clusters
TRsess=foreach (i = 10000:16060, .packages='survival', .combine='c') %dopar% {#getting threshold (via max R2):
  t1=as.numeric(substr(i,2,3))
  t2=as.numeric(substr(i,4,5))
  ifelse(t1+t2<59, {
    tr1=t1
    tr2=t1+t2
    de$time1=pmin(de$time,tr1)#spline functions for coefficients the slopes in each segment
    de$time2=pmax(tr1,pmin(de$time,tr2))
    de$time3=pmax(tr2,de$time)
    CLC_sess_=clogit(Va ~ Pl:time1+Ll:time1+Pl:time2+Ll:time2+Pl:time3+Ll:time3
                            +Pl_AUS+Ll_AUS+Pl_AUT+Ll_AUT+Pl_CAN+Ll_CAN+Pl_DEU+Ll_DEU+Pl_DNK+Ll_DNK+Pl_ESP+Ll_ESP+Pl_FIN+Ll_FIN+Pl_GBR+Ll_GBR+Pl_GRC+Ll_GRC+Pl_IRL+Ll_IRL+Pl_ISL+Ll_ISL+Pl_ISR+Ll_ISR+Pl_ITA+Ll_ITA+Pl_NLD+Ll_NLD+Pl_NOR+Ll_NOR+Pl_NZL+Ll_NZL+Pl_PRT+Ll_PRT+Pl_SWE+Ll_SWE
                            +strata(Esalt), data=de, method="efron")#regression
    CLC_sess_$loglik[2]
  } , NA)
}
stopCluster(cl)#stop parallel computing clusters
toc()
```

```
## 1203.204 sec elapsed
```

```r
tr=which.max(TRsess)
t1=as.numeric(substr((10000+(tr-1)),2,3))
t2=as.numeric(substr((10000+(tr-1)),4,5))
tr1=t1#9
tr2=t1+t2#37
tr1CLC_sess=tr1#
tr2CLC_sess=tr2#
tr1CLC_sess
```

```
## [1] 9
```

```r
tr2CLC_sess
```

```
## [1] 37
```

```r
#ALTERNATIVE SPECIFICATIONS SPLINE FOUR PERIODS:----

#optimization:
tic()
cores=detectCores()
cl=makeCluster(cores[1])
registerDoParallel(cl)#set parallel computing clusters
TRsesss=foreach (i = 1000000:1606060, .packages='survival', .combine='c') %dopar% {#getting threshold (via max R2):
  t1=as.numeric(substr(i,2,3))
  t2=as.numeric(substr(i,4,5))
  t3=as.numeric(substr(i,6,7))
  ifelse(t1+t2+t3<59, {
    tr1=t1
    tr2=t1+t2
    tr3=t1+t2+t3
    de$time1=pmin(de$time,tr1)#spline functions for coefficients the slopes in each segment
    de$time2=pmax(tr1,pmin(de$time,tr2))
    de$time3=pmax(tr2,pmin(de$time,tr3))
    de$time4=pmax(tr3,de$time)
    CLC_sesss_=clogit(Va ~ Pl:time1+Ll:time1+Pl:time2+Ll:time2+Pl:time3+Ll:time3+Pl:time4+Ll:time4
                         +Pl_AUS+Ll_AUS+Pl_AUT+Ll_AUT+Pl_CAN+Ll_CAN+Pl_DEU+Ll_DEU+Pl_DNK+Ll_DNK+Pl_ESP+Ll_ESP+Pl_FIN+Ll_FIN+Pl_GBR+Ll_GBR+Pl_GRC+Ll_GRC+Pl_IRL+Ll_IRL+Pl_ISL+Ll_ISL+Pl_ISR+Ll_ISR+Pl_ITA+Ll_ITA+Pl_NLD+Ll_NLD+Pl_NOR+Ll_NOR+Pl_NZL+Ll_NZL+Pl_PRT+Ll_PRT+Pl_SWE+Ll_SWE
                         +strata(Esalt), data=de, method="efron")#regression
    CLC_sesss_$loglik[2]
  } , NA)
}
stopCluster(cl)#stop parallel computing clusters
toc()
```

```
## 22591.699 sec elapsed
```

```r
tr=which.max(TRsesss)#
t1=as.numeric(substr((1000000+(tr-1)),2,3))
t2=as.numeric(substr((1000000+(tr-1)),4,5))
t3=as.numeric(substr((1000000+(tr-1)),6,7))
tr1=t1#13
tr2=t1+t2#
tr3=t1+t2+t3#37
tr1CLC_sesss=tr1#
tr2CLC_sesss=tr2#
tr3CLC_sesss=tr3#
tr1CLC_sesss
```

```
## [1] 13
```

```r
tr2CLC_sesss
```

```
## [1] 14
```

```r
tr3CLC_sesss
```

```
## [1] 37
```

```r
#ALTERNATIVE SPECIFICATIONS SPLINE FIVE PERIODS:----

#optimization:
tic()
cores=detectCores()
cl=makeCluster(cores[1])
registerDoParallel(cl)#set parallel computing clusters
TRsessss=foreach (i = 100000000:160606060, .packages='survival', .combine='c') %dopar% {#getting threshold (via max R2):
  t1=as.numeric(substr(i,2,3))
  t2=as.numeric(substr(i,4,5))
  t3=as.numeric(substr(i,6,7))
  t4=as.numeric(substr(i,8,9))
  ifelse(t1+t2+t3+t4<59, {
    tr1=t1
    tr2=t1+t2
    tr3=t1+t2+t3
    tr4=t1+t2+t3+t4
    de$time1=pmin(de$time,tr1)#spline functions for coefficients the slopes in each segment
    de$time2=pmax(tr1,pmin(de$time,tr2))
    de$time3=pmax(tr2,pmin(de$time,tr3))
    de$time4=pmax(tr3,pmin(de$time,tr4))
    de$time5=pmax(tr4,de$time)
    CLC_sessss_=clogit(Va ~ Pl:time1+Ll:time1+Pl:time2+Ll:time2+Pl:time3+Ll:time3+Pl:time4+Ll:time4+Pl:time5+Ll:time5
                         +Pl_AUS+Ll_AUS+Pl_AUT+Ll_AUT+Pl_CAN+Ll_CAN+Pl_DEU+Ll_DEU+Pl_DNK+Ll_DNK+Pl_ESP+Ll_ESP+Pl_FIN+Ll_FIN+Pl_GBR+Ll_GBR+Pl_GRC+Ll_GRC+Pl_IRL+Ll_IRL+Pl_ISL+Ll_ISL+Pl_ISR+Ll_ISR+Pl_ITA+Ll_ITA+Pl_NLD+Ll_NLD+Pl_NOR+Ll_NOR+Pl_NZL+Ll_NZL+Pl_PRT+Ll_PRT+Pl_SWE+Ll_SWE
                         +strata(Esalt), data=de, method="efron")#regression
    CLC_sessss_$loglik[2]
  } , NA)
}
stopCluster(cl)#stop parallel computing clusters
toc()
```

```
## 1175808.097 sec elapsed
```

```r
tr=which.max(TRsessss)#
t1=as.numeric(substr((100000000+(tr-1)),2,3))
t2=as.numeric(substr((100000000+(tr-1)),4,5))
t3=as.numeric(substr((100000000+(tr-1)),6,7))
t4=as.numeric(substr((100000000+(tr-1)),8,9))
tr1=t1#8
tr2=t1+t2#13
tr3=t1+t2+t3#14
tr4=t1+t2+t3+t4#37
tr1CLC_sessss=tr1#
tr2CLC_sessss=tr2#
tr3CLC_sessss=tr3#
tr4CLC_sessss=tr4#
tr1CLC_sessss
```

```
## [1] 8
```

```r
tr2CLC_sessss
```

```
## [1] 13
```

```r
tr3CLC_sessss
```

```
## [1] 14
```

```r
tr4CLC_sessss
```

```
## [1] 37
```

